# Author Omanjana Goswami for UCS analysis titled Less Fertilizer, Better Outcomes: USDA Conservation Programs Benefit Both Farmers and the Planet

# This research was performed in 2025, for report release in Feb 2026

#Data sources:

#Fertilizer use for US: Downloaded nitrogen fertilizer used in the US data from IFASTAT at https://www.ifastat.org/databases/graph/1_2 on August 11, 2025
#Fertilizer use for Illinois, Iowa, and Minnesota: Downloaded fertilizer used on corn (2021) and soybean (2023) data from USDA NASS Quickstats query tool on August 11, 2025
#Link for corn query: https://quickstats.nass.usda.gov/results/50DFA24D-16D1-3DF6-A345-82E67824EF9A#5812F8D6-EF01-3E0D-AAFA-AAE94B35DC47
#Link for soybean query: https://quickstats.nass.usda.gov/results/6ABC76A8-1866-3647-B779-3355AA65908E#CDFBB714-12FD-32BC-A335-64251BB38DAD

###########################################################################

#install packages
library(tidyverse)
library(dplyr)


#Set the working directory
getwd()
setwd("C:/Users/OGoswami/R files")

#import data into dataframe
#read csv in the folder to Global Environment
ndata <- read.csv("NitrogenApplied_ByState_OG.csv", header = TRUE, sep = ",")

#Confirm data in rows and columns

print(ndata)
colnames(ndata)
nrow(ndata)
ndata$NitrogenApplied <- as.numeric(gsub(",", "", ndata$NitrogenApplied))

# Convert nitrogen fertilizer applied on corn and soybean from pounds to metric ton

ndata %>%
  mutate(nitrogen_metrictons = NitrogenApplied * 0.00045359237)

# Calculate total nitrogen fertilizer applied by crop
Total_by_crop <- ndata %>%
  group_by(Crop) %>%
  summarise(Total_Nitrogen = sum(NitrogenApplied))

# Calculate total nitrogen fertilizer applied by each state
Total_by_state <- ndata %>%
  group_by(State) %>%
  summarise(Total_Nitrogen = sum(NitrogenApplied))

# Print results of aggregation where nitrogen fertilizer is in pounds 
print(Total_by_crop)
print(Total_by_state)

#convert units of nitrogen applied in each state from pounds to metric tons using the conversion factor of 0.00045359237

ndata %>%
  group_by(State) %>%
  summarise(Total_Nitrogen = sum(as.numeric(gsub(",", "", NitrogenApplied)),
      na.rm = TRUE),
    .groups = "drop") %>%
  mutate(nitrogen_metrictons = Total_Nitrogen * 0.00045359237) %>%
  mutate(Excess50 = nitrogen_metrictons* 0.5) %>%
  mutate(Excess30 = nitrogen_metrictons* 0.3) %>%
         select(State, nitrogen_metrictons, Excess50, Excess30)
  
##########################################################################
# This section will calculate the excess fertilizer use nationally

#According to 2023 data from IFASTAT 11620.1 thousand tonnes (=11.620100 million metric ton) of N fertilizer was used  

# To calculate the quantity of nitrogen fertilizer that is excess assuming 50 percent NUE

N_fert_national <- 11620100
Excess50_fert_national_metrictons<- N_fert_national* 0.5

print (paste(N_fert_national, "metric ton"))
print (paste(Excess50_fert_national_metrictons, "metric ton"))

# To calculate the quantity of nitrogen fertilizer that is excess assuming 70 percent NUE
# NUE of 70 percent means that of the quantity applied 70 percent is used by the plant and the remaining 30 percent is excess
Excess30_fert_national_metrictons<- N_fert_national* 0.3

print (paste(N_fert_national, "metric ton"))
print (paste(Excess30_fert_national_metrictons, "metric ton"))


################################################################################

# This section of the code will now calculate the heat trapping gas emissions for US and states from fertilizer use using the IPCC equations
# The definition, source, and unit of each term is available in the accompanying report (Appendix A)
# We focus on emissions of Nitrous oxide (N2O) and carbon dioxide (CO2) from breakdown of fertilizer and liming
# All computed values are in units of metric ton CO2 equivalent

# Equation 1: EdirectN2O = FN * EF1 * (44/28) * 273
# Equation 2: EGASF= FN * FracGASF * EF4 *  (44/28)  * 273					     	
# Equation 3: ELEACH= FN * FracLEACH * EF5 *  (44/28)  * 273				
# Equation 4  Edecomposition  = FN * EFdecomposition				
# Equation 5  Elimestone = FN* LT * 0.12 *  44/12

# The definition of the terms used in the above equations are as follows: 
# FN = Amount of applied nitrogen fertilizer in Mt N
# EF1 = 0.01 from IPCC Chapter 11 
# FracGASF = Fraction of nitrogen fertilizer that volatilizes as ammonia and NOx
# EF4 =  Emission factor for N2O from atmospheric deposition 
# FracLEACH = Fraction of nitrogen fertilizer that leaches as nitrate
# EF5 = Emission factor for N2O from leached nitrate
# EFdecomposition = conversion factor of embedded carbon in urea which is 1.57 (urea is the most common fertilizer used in the US)
# LT = demanded limestone to neutralize the soil per ton of applied nitrogen


#read csv with data in the folder

Equation_Data <- read.csv("FertUse_GHGcalc_OG.csv", header = TRUE, sep = ",")
print(Equation_Data)

# Assigning values to the variables
EF1 <- 0.01
FracGASF <- 0.11
EF4 <- 0.01
FracLEACH <- 0.24
EF5 <- 0.011
EFdecomposition <- 1.57
LT <- 6.5

#Creating a new dataframe and adding columns with calculations

Equation_Data <- Equation_Data %>%
  mutate (FN = as.numeric(Fertilizer.Applied.in.Metrictons))%>%
  mutate (EdirectN2O_metrictons = FN * EF1 * (44/28) * 273)%>%
  mutate (EGASF_metrictons = FN * FracGASF * EF4 *  (44/28)  * 273)%>%
  mutate (ELEACH_metrictons = FN * FracLEACH * EF5 *  (44/28)  * 273)%>%
  mutate (Edecomposition_metrictons  = FN * EFdecomposition)%>%
  mutate (Elimestone_metrictons = FN * LT * 0.12 *  (44/12)) %>%
  mutate (Total = EdirectN2O_metrictons + EGASF_metrictons + ELEACH_metrictons + Edecomposition_metrictons + Elimestone_metrictons)

# To save the calculations in tabular format in a new CSV 
write.csv(Equation_Data, "GHG_eqtn_values_state_national_OG.csv")

################################################################################

# This section of the code will now calculate heat trapping emissions for US and states from excess fertilizer use (assuming 50% NUE)
# The same IPCC equations are being used here
# The definition, source, and unit of each term is available in the accompanying report (Appendix A)
# We focus on emissions of Nitrous oxide (N2O) and carbon dioxide (CO2)
# All computed values are in units of metric ton CO2 equivalent

#read csv containing data in the folder to Global Environment

Equation_Data_for_excess <- read.csv("Excess50_FertUse_GHGcalc_OG.csv", header = TRUE, sep = ",")
print(Equation_Data_for_excess)

EF1 <- 0.01
FracGASF <- 0.11
EF4 <- 0.01
FracLEACH <- 0.24
EF5 <- 0.011
EFdecomposition <- 1.57
LT <- 6.5

#Creating a new dataframe and adding columns with calculations

Equation_Data_for_excess <- Equation_Data_for_excess %>%
  mutate (FN = as.numeric(Excess.fertilizer.Applied.in.Metrictons))%>%
  mutate (EdirectN2O = FN * EF1 * (44/28) * 273)%>%
  mutate (EGASF = FN * FracGASF * EF4 *  (44/28)  * 273)%>%
  mutate (ELEACH = FN * FracLEACH * EF5 *  (44/28)  * 273)		%>%
  mutate (Edecomposition  = FN * EFdecomposition)%>%
  mutate (Elimestone = FN * LT * 0.12 *  (44/12)) %>%
  mutate (Total = EdirectN2O + EGASF + ELEACH + Edecomposition + Elimestone)

# To save the calculations in tabular format in a new CSV
write.csv(Equation_Data_for_excess, "Excess50_GHG_eqtn_values_state_national_OG.csv")

################################################################################### 

# This section of the code will now calculate heat trapping emissions for US and states from excess fertilizer use (assuming 70% NUE)
# The same IPCC equations are being used here
# The definition, source, and unit of each term is available in the accompanying report (Appendix A)
# All computed values are in units of metric ton CO2 equivalent

#read csv containing data in the folder

Equation_Data_for_excess <- read.csv("Excess30_FertUse_GHGcalc_OG.csv", header = TRUE, sep = ",")
print(Equation_Data_for_excess)

EF1 <- 0.01
FracGASF <- 0.11
EF4 <- 0.01
FracLEACH <- 0.24
EF5 <- 0.011
EFdecomposition <- 1.57
LT <- 6.5

#Creating a new dataframe and adding columns with calculations

Equation_Data_for_excess <- Equation_Data_for_excess %>%
  mutate (FN = as.numeric(Excess.fertilizer.Applied.in.Metrictons))%>%
  mutate (EdirectN2O = FN * EF1 * (44/28) * 273)%>%
  mutate (EGASF = FN * FracGASF * EF4 *  (44/28)  * 273)%>%
  mutate (ELEACH = FN * FracLEACH * EF5 *  (44/28)  * 273)		%>%
  mutate (Edecomposition  = FN * EFdecomposition)%>%
  mutate (Elimestone = FN * LT * 0.12 *  (44/12))%>%
  mutate (Total = EdirectN2O + EGASF + ELEACH + Edecomposition + Elimestone)
  
# To save the calculations in tabular format in a new CSV
write.csv(Equation_Data_for_excess, "Excess30_GHG_eqtn_values_state_national_OG.csv")

#####END
